clear
set more off

local base = 1

insheet using qtr_means_heroin_pain.csv, comma clear
sort year quarter
keep if year>=2004
gen trend=4*(year-2004)+quarter
xi i.quarter
global mint = 19
global maxt = 36


tempfile main bootsave 

qui save `main' , replace 

if `base' == 1 {

	postfile bskeep trend sseh sseo sseol Fh Fp Fpl using quadratic_spline_results_pain, replace


	forvalues b = $mint/$maxt { 
		display `b'
		use `main', replace 
		gen z=trend-`b'
		gen post=trend>=`b'
		gen zpost1=post*z
		gen zpost2=zpost1*zpost1
		gen zpre1=(1-post)*z
		gen zpre2=zpre1*zpre1
		reg heroin_30 zpre1 zpre2 zpost1 zpost2
		local sseh=e(rss)
		test (zpre1=zpost1)(zpre2=zpost2)
		local Fh = r(F)
		reg pain_30 zpre1 zpre2 zpost1 zpost2
		local sseo=e(rss)
		test (zpre1=zpost1)(zpre2=zpost2)
		local Fp = r(F)
		reg pain_30 zpre1 zpost1 
		local sseol=e(rss)
		test (zpre1=zpost1)
		local Fpl = r(F)

		post bskeep (`b') (`sseh') (`sseo') (`sseol') (`Fh') (`Fp') (`Fpl')
	}
	postclose bskeep
	clear
	use quadratic_spline_results_pain
	list
	clear


	use `main', replace
	sort trend
	merge 1:1 trend using quadratic_spline_results_pain
	egen minsseh=min(sseh)
	egen minsseo=min(sseo)
	egen minsseol=min(sseol)

	gen mh=minsseh==sseh
	gen mo=minsseo==sseo
	gen mol=minsseol==sseol

	egen th=max(mh*trend)
	egen to=max(mo*trend)
	egen tol=max(mol*trend)


	gen z=trend-th
	gen post=trend>=th
	gen zpost1=post*z
	gen zpost2=zpost1*zpost1
	gen zpre1=(1-post)*z
	gen zpre2=zpre1*zpre1
	reg heroin_30 zpre1 zpre2 zpost1 zpost2
	test (zpre1=zpost1) (zpre2=zpost2)
	predict herion_30p
	drop z post zpost* zpre*

	gen z=trend-to
	gen post=trend>=to
	gen zpost1=post*z
	gen zpost2=zpost1*zpost1
	gen zpre1=(1-post)*z
	gen zpre2=zpre1*zpre1
	reg pain_30 zpre1 zpre2 zpost1 zpost2
	test (zpre1=zpost1) (zpre2=zpost2)
	predict pain_30p
	drop z post zpost* zpre*

	gen z=trend-tol
	gen post=trend>=tol
	gen zpost1=post*z
	gen zpre1=(1-post)*z
	reg pain_30 zpre1 zpost1 
	test (zpre1=zpost1)
	predict pain_30pl
	drop z post zpost* zpre*
	summ tol
	local breaks = r(mean)
	twoway (line pain_30 trend)(line pain_30pl trend), xline(`breaks')
	
	keep year quarter trend heroin_30 herion_30p sseh th pain_30 pain_30p sseo Fh Fp
	order year quarter trend heroin_30 herion_30p sseh th pain_30 pain_30p sseo Fh Fp
	outsheet using quadratic_spline_results_pain.csv, comma replace
}

